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Abstract 


Atomistic simulations of intergranular fracture have indicated that grain-scale crack 
growth in polycrystalline metals can be direction dependent. At these material length 
scales, the atomic environment greatly influences the nature of intergranular crack 
propagation, through either brittle or ductile mechanisms, that are a function of adjacent 
grain orientation and direction of crack propagation. Methods have been developed to 
obtain cohesive zone models (CZM) directly from molecular dynamics simulations. 
These CZMs may be incorporated into decohesion finite element formulations to 
simulate fracture at larger length scales. A new directional decohesion element is 
presented that calculates the direction of Mode I opening and incorporates a material 
criterion for dislocation emission based on the local crystallographic environment to 
automatically select the CZM that best represents crack growth. The simulation of 
fracture in 2-D and 3-D aluminum poly crystals is used to illustrate the effect of 
parameterized CZMs and the effectiveness of directional decohesion finite elements. 
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1. Introduction 


At nanometer length scales, the local atomic environment greatly influences the nature 
of intergranular crack propagation. Crack tip process zones can be characterized by brittle 
or ductile mechanisms that are initiated as a function of adjacent grain orientation and 
direction of crack propagation. Experimental evidence for the crystallographic 
dependence of crack tip processes date back to the work of Bilby and Bullough (1954). 
The atomistic response at the crack front can be obtained through molecular dynamics 
(MD) simulation and cast into an aggregate traction-displacement relationship and 
subsequently used as part of a sequential multiscale analysis (Yamakov et al., 2006; 
Saether, 2008). 

Sequential multiscale modeling typically involves some form of averaging of physical 
parameters that can serve as initial conditions or provide material parameters to another 
model which is analyzed separately. A desirable aspect of sequential methods is that 
length and time scales between independent material models do not have to be coupled. 
The averaging or homogenization of information across length scales that is inherent to 
sequential multiscale methods is depicted in Figure 1 where a notional coupling is shown 
across domains representing characteristic features at the sub-atomic through structural 
length scales. 

Cohesive zone models (CZM) are based on traction-displacement relationships and are 
a viable tool for sequential multiscale modeling (Saether, 2008). Within a sequential 
analysis framework, CZMs based on averaged atomistic deformation processes for 
intergranular and transgranular fracture can be used to simulate material behavior at 
larger length scales. 

The notional flow of sequential coupling that utilizes CZMs to carry information of 
microscopic failure mechanisms to predict damage progression at larger length scales is 
depicted in Figure 2. Hitherto, the properties of the CZM have been determined 
empirically or heuristically, and mixed-mode fracture has been limited to empirical 
relationships that have been applied at both microscopic and macroscopic length scales. 
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In Figure 2, CZMs provide the critical transition between inherently atomistic and 
inherently continuum representations. 
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Figure 1. Hierarchy of models over length scales from quantum mechanical 
through homogeneous continuum. (Adapted from Oden et al., 2006) 
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Figure 2. Multiscale analysis with cohesive zone models. 
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Characterization of intergranular fracture for different grain boundary (GB) 
orientations is currently very limited. Yamakov et al. (2006) determined CZM parameters 
for intergranular fracture of a symmetric tilt Z99 GB, which, depending on the direction 
of propagation, exhibited either ductile or brittle fracture characteristics in Mode I 
opening. With the lack of predictions of CZM parameters for other GB orientations, 
scaling the results obtained by Yamakov et al. (2006) for fracture along a single GB is 
used in the current study to represent other GB orientations. 

CZMs may be incorporated into decohesion finite element formulations to simulate 
fracture at continuum length scales. The CZMs provide constitutive relations that govern 
the separation of material planes due to relative displacement measures that correspond to 
the three fundamental modes of fracture. A new feature to decohesion finite element 
formulations is presented that allows the element to automatically apply appropriate 
brittle or ductile CZMs for crack growth depending on the direction of fracture 
propagation. This new feature thus modifies the element behavior as the finite element 
analysis proceeds. 

Six-node and 12-node directional decohesion finite elements are formulated that 
incorporate Rice’s criterion (Rice, 1992) for dislocation emission to automatically select 
the appropriate CZM to simulate brittle or ductile crack propagation characteristics. For 
intergranular fracture, element input requires the crystallographic orientation of the 
adjacent grains in the microstructure and additional material parameters including 
surface, GB, stable and unstable stacking-fault, and unstable twinning energy in order to 
predict brittle or ductile Mode I fracture. 

Section 2 presents a general discussion of cohesive zone models, and Section 3 
presents an overview of decohesion element formulations. Specific formulations of 
directional decohesion elements are detailed in Section 3.1. Finally, Section 4 presents 
illustrative simulations of intergranular fracture in 2-D and 3-D aluminum polycrystal 
models. A 2-D model is used to demonstrate the effect of CZM parameterization, and two 
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3-D models are used to assess the unique modeling capability of directional decohesion 
finite elements. 


2. Cohesive Zone Models 

Cohesive zone models were originally developed to represent complicated nonlinear 
fracture processes in ductile and quasi-brittle materials (Dugdale, 1960; Barenblatt, 
1962). CZMs were later developed to describe general adhesion and frictional slip along 
an interface (Maugis, 1992; Kem et al., 1998). The CZM approach is formulated from a 
constitutive relationship based on applied tractions and relative displacements to 
represent separation in various fracture modes of two initially coincident and bonded 
internal surfaces. The relative displacements associated with the creation of a new 
fracture surface for the three fundamental fracture modes are depicted in Figure 3. 



Figure 3. Fundamental fracture modes in solids 


In a CZM-based approach, the area under the CZM traction-displacement curve 
represents the work of separation to open a crack in a particular fracture mode, i.e. the 
fracture toughness, and is given by Tvergaard and Hutchinson (1992) as 

A f 

G c = JrdA (1) 

0 
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where G c is the work of separation, r is the applied traction, A is a general form of the 
displacement, and A f is the critical displacement at which complete separation has 
occurred and the tractions are zero. 

In general, the traction-displacement relationships t(a) are obtained through 
differentiation of a potential (/)=(f){A), which represents the free energy of decohesion. 
The selection of a potential function is typically based on recovering the assumed 
traction-displacement relationship, and particular forms are generally selected for 
analytical convenience. In practice, various forms have been used as shown in Table 1. 
The existence of a work potential yields the work of separation as independent of the 
loading history and of the specific shape of the traction-displacement law. 

Cohesive properties along an interface have typically been approximated using 
empirical data to define the CZMs (Tvergaard and Hutchinson, 1992; Costanzo and 
Allen, 1995; Camacho and Ortiz, 1996; Klein and Gao, 1998; Zavattieri et al., 2001; 
Zavattieri and Espinosa, 2003; Turon et al., 2004). These models are frequently used in 
conjunction with the finite element method (FEM) to study fracture at macroscopic 
length scales in a wide variety of materials. 

The shape of the CZM law represents the behavior of the material near the crack 
tip under load. Various attempts have been made to determine the shape based on 
fundamental bonding characteristics in metals (Rose et al., 1983; Nguyen and Ortiz, 
2002). The most commonly assumed forms of the traction-displacement law have been 
expressed as exponential, bilinear, and trapezoidal functions. A general review of 
various forms of CZMs is presented by Chandra et al. (2002) and summarized in Table 1. 
The table illustrates the basic form of the CZMs, their key parameters, and important 
features. 
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Table 1. Various cohesive zone models and their parameters. (From Chandra et al., 2002) 


Author 
(year) 
Barenblatt 
(1959, 1962) 


Dugdale 

(1960) 


Needieman 

(1987) 


Rice and 

Wang 

(1989) 


Needieman 

(1990a) 


Needieman 

(1990b) 


Tvcrgaard 

(1990) 


Tvcrgaard 
and Hutch- 
inson (1992) 


Xu and 

Needieman 

(1993) 


Proposed model 



Model parameters 

Problem 

solved 

Model constants 

Comments 

tr i^ +et ' g i(Qdf 

A Jo s 

— (ductile) 

- (brittle) 

Perfectly 

brittle 

materials 


The first to 
propose the 
cohesive zone 
concept 

t = r 0 + 7i 

7o is work of separation 
for brittle materials 
T { is work of plastic deformation 




~ = 2 sin 2 

For small value of T/a y 
f = L23^ 2 

Yielding of 
thin ideal 
elastic-plastic 
steel sheets 
containing slits 

Plastic zone 
ranges from 
0.042 to 0.448 
(in.) 

Cohesive stress 
equated to yield 
stress of material 




4> sep is work of separation 
S are normalizing parameters 
Cmax is cohesive strength 


Particle-matrix 

8=\0~ 9 to 

Phenomenological 

decohesion 

10“ 8 m cohesive 

model; predicts 


energy 1-10 J/m 2 
- 1000 

normal separation 


1400 MPa 



ff y = 350-450 MPa 




Model based on 
atomic fit of the 
type 0+x)T 
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E 0 is initial Young’s modulus 
h is normalizing parameter 
0Wx is maximum stress 
a is constant^ = 2y) 


Solute 

segregation 


Ascending part is 
equated to E 0 
considers normal 
separation and 
ignores shear 
separation 


Polynomial 


Exponential fit 
for T n . Linear 



(psep is work of separation 
8 are normalizing parameters 
cr max is cohesive strength 


Particle-matrix 8 — 10 9 to 
decohesion 1CT 8 m 


Predicts normal 
separation 


Exponential fit for normal T n 
Trigonometric fit for shear r, 




<fr n . <j> t arc work of normal 

Deeohesion of 


Periodic shear 

and shear separation 

interface under 

2 x 10 ,0 to 

traction to model 

d„, d, arc critical displacements 

hydrostatic 

2 x 10-’ m 

Pieriels shear 

<r UM is cohesive strength 

tension 

■Ac/*. = 
0.57-2.59 
<w/*e — 2, 3 

stress due to slip 


fi„, A are critical displacements 
is cohesive strength 

Interfaces of 
whisker rein- 
forced metal 
mutrix com- 
posites 

1 x JO' 9 m 
E - 60 GPa 
(Young’s mod) 
o y /E 0.005 
<W/*y ~ 5-9 

Quadratic model 

is work of separation 

Crack growth 

r„/r 0 o-io 

Claims shape of 

<\ is critical displacement 

in etaslo-plas- 

(/.= Plwk.), 

separation law are 

Cn»x is peak normal traction 

tic material, 

W * i 

relatively unim- 

/interface strength d|* 

Sj arc factors governing shape 

peeling of ad- 
hesive joints 

GmxxfOy 0— I4d|, 

-0.15. 0.5 
Oy/E- 1/300 

portant 


Exponential fit for both T n and Tf 


i />„ is work of normal separation 
4>< is work of shear separation 
«>„. 8, are critical displacements 
<r lna( is cohesive strength 


Particle-matrix - 8, — 
deeohesion 2 x 10" 10 to 

2 x 10 9 m 


Predicts shear and 
normal separation 



Despite all of the forms of CZMs that have been proposed, a commonly used 
mathematical form not shown in Table 1, is a bilinear constitutive relation. The bilinear 
form is often chosen because of its mathematical simplicity and its suitability for 
representing brittle and ductile fracture in metals (Yamakov et al., 2006), and purely 
brittle fracture in polymeric and ceramic composite materials (Camacho and Ortiz, 1996). 
Bilinear types of cohesive zone models have been used in several recent finite element 
simulations of brittle fracture during multi-axial dynamic loading of ceramic 
microstructures (Zavattieri et al., 2001, 2003). This form of a CZM is shown in Figure 4. 



In Figure 4, zl represents relative discontinuous displacement jumps between the 
upper and lower nodes, A = A top - A 1 ’ 0 '. A° indicates the displacement jump corresponding 
to the maximum traction and indicates the onset of damage. Final failure of the cohesive 
zone is assumed after the relative displacement jump yields zero traction. This maximum 
relative displacement jump is indicated as A. 



An interpolation procedure using single-mode CZMs has been developed by Turon et 
al. (2004) to generate an effective CZM for mixed-mode applications. This interpolation 
or coupling scheme is depicted in Figure 5. 


T 



In Figure 5, the relative displacement jumps for single -mode CZMs are denoted by 
S = 5 t0 p - S bot while A refers to the interpolated displacement jumps. The interpolation of 
“normal” and “shear” CZMs to obtain a single effective CZM is made to represent mixed 
normal and tangential fracture modes. The Mode II and III displacements, Si and Sw, are 
combined by the root sum of the squares to make a combined “shear” displacement, 

8 S = V < 8 . i ) 2 +( 5 ...) 2 ( 2 ) 

and the coupled displacement across the interface is defined by 

A = ,/(S,} 2 +(8 s ) 2 (3) 
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where ( ) is the MacAuley bracket function, (x) = A (x + |x|). Mode mixity parameters 
are defined by 


P = 


^_s 

8s+(8 I > 


B = 


P 2 

1 + 2p 2 - 2(3 


(4) 


where J3 = 0 for pure Mode I loading, and J3 = 1 for pure shear loading. The fracture 
onset criterion is defined by 


A 0 



B n 


(5) 


and the criterion for complete fracture representing crack propagation is derived from the 
B-K critical energy release rate expression (Benzeggagh and Kenane, 1996) given by 


Gc= G lc + (G shear - G Ic ) 


G 


shear 


v Gtot J 


( 6 ) 


where G to/ = G\ + G s hear, G s h ea r = Gw + Gm, and r/ is an empirical factor used to correlate 
with experimental data. This leads to an expression for the final opening displacement 
jump corresponding to crack propagation given by 

j 8f 5f + (8g5g - S^Sf ) B n 

A 0 

In the previous expressions, df , df, , and df n are the individual nodal displacement 
jumps corresponding to opening in Mode I, II and III, respectively, and are user- 
specified. The S° , djj , and <?m are the individual nodal displacements corresponding to 
the onset of damage and are computed from the user-specified initial stiffness values, , 
£jj , and k ' m , and peak tractions, t ° , t\ j , and t ° m , of the single-mode CZMs. The 
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interpolated parameters, A° and A, define the mixed-mode CZM used in decohesion 
finite elements for the continuum simulation of fracture. 

2.1 CZM Parameters Dependent on Direction of Crack Growth 

Important damage and failure behavior in polycrystalline metals can be simulated 
including intergranular fracture along grain boundaries in metals, transgranular fracture 
along slip planes, and fracture through interfaces between grains and small second phase 
particles. At atomistic length scales, it has been shown that for a common 299 symmetric 
tilt GB, Mode I fracture can propagate by distinct ductile or brittle modes depending on 
the direction of propagation (Yamakov et al., 2006). A synopsis of the distinct directional 
difference in crack propagation behavior and the atomistic determination of 
corresponding CZM parameters as presented by Yamakov et al. (2006) is shown in 
Figure 6. The system was subjected to a hydrostatic tensile load in order to avoid shear 
stresses away from the crack that could induce undesired dislocation formation. To 
ensure crack nucleation and growth at the nanometer scale, the prestress in the range of 
3.75 to 4.25 GPa was considered. 

Figure 6 shows characteristics of the propagating internal crack at different times 
during a MD simulation. At t = 12 ps, the center crack shown at the top of the figure has 
well-developed differences in the crack tip behavior at the two ends. The left crack tip, 
highlighted at t = 68 ps, shows a large degree of blunting due to several formations that 
are identified as the development of twinning (©), cores of partial or twinning 
dislocations ((D) and secondary slip (©). The right crack tip, shown at t = 55 ps, 
demonstrates cleavage associated with brittle fracture. The calculated CZM for ductile 
opening shows a peak traction of ~4.0 GPa and an expected large opening displacement 
for final fracture. For brittle fracture, the CZM is characterized by a peak traction of ~5.0 
GPa and a shorter softening region prior to complete opening in Mode I. The CZM 
curves are computed as moving averages of traction-displacement measurements 
obtained from MD simulation (See Yamakov et al., 2006). In the CZM curves presented 
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in Figure 6, the individual data points and a least-squares curve fit for a prestress of 4.25 
GPa are shown, together with curve fits for the other magnitudes of applied prestress. 



Ductile crack tip 


Brittle crack tip 



-20 -10 r , G 10 20 

x[aj 



Ductile CZM 


— 4.25 GPk 

— 4.25 GP& 

— 4,0 GPa 

5 ( il',i 


Brittle CZM 



Figure 6. Direction-dependant crack growth behavior incorporated into 
decohesion finite element formulations (From Yamakov et ah, 2006). 
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The brittle and ductile response of the crack tips shown in Figure 6 is a function of 
several parameters that include the atomic configuration (e.g. the orientation of slip 
planes on either side of the grain boundary), and the magnitude of applied external loads. 

Figure 7 shows a configuration of a crack plane and slip plane with the angle between 
the two specified by 0. The orientation of the Burgers vector, b, lying in the slip plane, 
showing the direction of slip is given by the angle ^ . For a slip direction aligned with the 
crack opening such that <j>= 0, the Rice criterion (Rice, 1992) for dislocation emission 
associated with Mode I opening along two propagation directions can be written as 

_ Yus _ (l±cos0)sin 2 0 a > 0 for brittle fracture ™ 

a 2y s -y GB 8 W ere a < 0 for ductile fracture 

where + cos# is associated with the -x direction and cos(/r — d) = -cos 6 is associated 
with the +x direction, and a is the numerical value of the criterion used to predict brittle 
or ductile fracture depending on its sign. For a 299 symmetric tilt GB in aluminum, 
Yamakov et al. (2006) identifies an unstable stacking-fault energy of y us = 0.168 J/m 2 , a 



Figure 7. Geometric configuration of a crack with an intersecting slip plane. 
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surface energy of the GB as y s = 0.952 ±0.010 J/m 2 , and an excess GB energy given by 
Ygb = 0.60 ± 0.05 J/m 2 . Substituting in the energy quantities and multiplying through by 
the denominator in the second term of Equation (8) (because only the sign of a is 
important), the resulting Rice criterion for this system can be expressed as 

a = 1 .0307 ± 0.01 - (l ± cos 0)sin 2 0 . (9) 


3. Decohesion Finite Elements 

Decohesion finite elements (also referred to in the literature as “interface elements”, 
“cohesive elements,” or “CZM elements”) were initially introduced by Hilleborg et al. 
(1976) to study fracture in concrete. Needleman (1987, 1990a, 1990b) later formulated a 
decohesion element to study dynamic fracture in isotropic solids. The formulations of the 
elements presented here are based on the work of Beer (1985). The basic topology of 
decohesion elements is formed by separable surfaces that can open according to 
constitutive relationships that relate tractions and displacements. The element 
formulations have found broad application for modeling brittle and ductile fracture 
whether or not crack paths are known a priori. 

When the location of fracture is not known a priori, decohesion elements can be 
placed between all continuum elements within a finite element mesh. While decohesion 
elements can add artificial compliance to the model (Diehl, 2008), more complex 
simulations using adaptive meshing can limit the number of decohesion elements by 
placing these elements only in regions of high stress (Zhang et al., 2007). Figure 8 depicts 
a micromechanical simulation of intergranular fracture with decohesion elements placed 
along all grain boundaries. 

The formulation of decohesion finite elements is based on the dimensionality, order of 
element interpolants, and the particular choice of CZMs that are incorporated into the 
formulation. Decohesion elements are formulated with two sets of initially coincident 
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nodes defining two superposed surfaces forming a cohesive interface. In Figure 9, ff 
represents the middle surface of the element, and £T and Q represent the top and bottom 
surfaces of the element after separation. The particular configuration depicted in Figure 9 
represents an 8-node linear decohesion element used to simulate the cohesive behavior 
between connected 8-node solid continuum elements. Shape functions appropriate to the 
order of assumed variation over the element domain govern the interpolation of quantities 
over the superposed surfaces. 




Figure 8. Embedding decohesion elements along GBs to study microstructural fracture. 



Figure 9. Decohesion element cohesive surfaces. 
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Fundamental to the response of decohesion elements is the form of the assumed 
cohesive zone models (CZMs) that govern the prediction of separation - fracture - of 
element surfaces (Rabinovitch, 2008). As discussed in Section 2, CZMs are defined by 
displacement jumps, denoted by A t , which measure the discontinuity of the i ,h 
displacement components between initially coincident top ( S + ) and bottom points on the 

surfaces ( S ~ ). Thus, the resulting finite element is initially of zero thickness across the 

assumed fracture plane (see Beer, 1985; Camanho et al., 2003; and Cook et al., 2001). 
Discretizing each surface into k finite element nodes and assuming that the displacement 
field varies over the surface defined by standard Lagrangian shape functions Nt yields the 
relationship 


[A,] 



- s( 


N k 8 lk 

^2 

> = < 

82 

- 5 - 

> = < 

N k S 2k 

^3 


83 

-83 


N k 8 3k 


( 10 ) 


where is the i th relative displacement at the k ,h node referenced to the middle surface of 
the element /2°. Transforming between local ((?/,e?,e,;) and global (x/.x^xj) coordinate 
systems using the transformation tensor Ty yields 


{u} = [ tKa} = [t][n]{8} = [bKs} (11) 


Relating the element tractions, r„ to the relative element displacements, u h is performed 
using the decohesion element constitutive operator, Dy, as 

{t} = [d]{u} (12) 


The governing equation can be obtained from the principle of virtual work 


Jd{u} T {x}dQ° - {f} T d{u} = 0 (13) 

Q° 
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which yields, for a geometrically linear problem, the relationship 


\[B] T {T}df2° -{f} T =0 (14) 


where [5] is constant. Substituting Equations (11) and (12) into Equation (14) yields the 
standard governing equation for linear elastostatics 


[K]{S}={f} (15) 

where 

[K]= |[B] T [D][B]dQ° (16) 

Q° 

The constitutive matrix, \D], is formulated based on the assumed traction-displacement 
relationship. 

Different features can be incorporated into the element formulation. For example, 
Turon et al. (2004) utilizes a tangent constitutive relation based on bilinear CZMs that 
incorporates a damage parameter as a state variable for mixed-mode failure and a 
variable term that accounts for interpenetration due to possible closing in Mode I fracture. 
Various additional details of decohesion element formulations can be found in Chen et 
al., 1999; Foulk et al., 2000; Alfano and Crisfield, 2001; Goyal et al., 2002; Camanho et 
al., 2003; Segurado and Llorca, 2004; Gustafson et al., 2008; Hamitouche et al., 2008; 
and Van den Bosch et al., 2007. However, most of the differences in these formulations 
pertain to specific definitions of CZM relationships incorporated into the decohesion 
element. 

An additional issue relevant to the development of CZMs regards the initial stiffness 
associated with the constitutive relationship for each fracture mode. For macroscopic 
analysis, this initial stiffness assumes the role of a penalty constraint. A formal penalty 
parameter would assume the highest value that would not cause ill-conditioning in the 
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global element stiffness matrix. However, in practice, this constraint is set to a high value 
to minimize the relative displacements between the top and bottom surfaces before the 
onset of softening, but does allow numerically small displacements to be calculated that 
are compared with the onset relative displacement at which softening begins. This onset 
displacement, denoted by zl 0 in Figure 4, is equal to the peak traction divided by the 
selected penalty stiffness and is, thus, a small quantity. However, a consequence of this is 
the introduction of a spurious compliance into the model that can alter results, 
particularly in models where decohesion elements are placed between many continuum 
elements. However, at nanoscopic length scales, this initial stiffness is typically 
calculated based on the elastic properties of the grain boundary or between lattice planes 
because the thickness of the cohesive surface can no longer be generally regarded as 
infinitesimal. In this case, the initial stiffness of the CZM no longer approximates a 
mathematical penalty parameter but is reduced to represent the actual stiffness of a finite 
thickness cohesive zone. Under arbitrarily small loads, the displacement field is not 
homogeneous, and the cohesive surface exhibits an elastic separation due to straining 
below the peak traction. To avoid adding spurious compliance to the model due to a 
numerically finite initial stiffness, either an initial thickness may be assigned to the 
decohesion elements or the elastic properties of adjacent continuum elements may be 
adjusted. 

Convergence difficulties have been encountered and have been related to numerical 
problems caused by abrupt changes in the slope of the CZM relationship - such as at the 
vertices of triangular and trapezoidal CZM relations - during loading and unloading 
cycles in which large variations can occur in computing the tangent stiffness matrix (Gao 
and Bower, 2004). Another effect that can hinder convergence involves the type of 
numerical quadrature used to evaluate the element integrals. It has been found that 
Gaussian quadrature tends to couple kinematic degrees of freedom across the element 
that can be seen in the eigenmodes of the element deformation states. Conversely, 
Newton-Cotes quadrature methods act to uncouple the eigenmodes and have been 
associated with an improvement in convergence behavior (Schellekens and de Borst, 
1993). Additional numerical aspects, such as bounds on element size, have been 
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presented (Allen and Searcy, 2000; Tomar et al., 2004; and Turon et al., 2007), and 
augmented solution schemes that include numerical relaxation to improve convergence 
characteristics have been advanced (Gao and Bower, 2004). 


3 . 1 Directional Decohesion Element Formulation 

The directional dependence of Mode I fracture presented in Section 2.1 motivated the 
development of a directional decohesion finite element formulation for simulating 
intergranular crack propagation. The unique feature of this element is the adaptive 
application of different CZMs depending on the grain orientation and the internal 
computation of the direction of Mode I propagation. This element incorporates the 
direction-dependant CZMs for Mode I opening and utilizes a single Mode II and Mode 
III CZM for computing an effective mixed-mode CZM as explained in Section 2.0. The 
appropriate brittle or ductile Mode I CZM is selected based on Rice’s criterion that 
accounts for the orientation of the surrounding slip systems in the direction of 
propagation (see Figure 6). Both 1-D and 2-D element configurations are presented 
which incorporate Rice’s criterion (Rice, 1992) to apply the appropriate CZM during an 
incremental-iterative solution procedure. The directional decohesion element 
formulations have been encoded into user-defined element (UEL) subroutines to enable 
simulations to be performed using the commercial finite element code 
ABAQUS® 1 /Standard (2008). 

Until the final effective opening displacement has been reached and the decohesion 
element stiffness terms are zero, calculations are performed using the relative u, v, and w 
nodal displacements to approximate the displacement gradients over the element surface 
and determine the direction of propagation. For cases in which various opening 
displacements are equal or zero, a default direction is assumed. A propagation vector, d , 
is used to specify the direction of Mode I fracture as is required in the application of 
Rice’s criterion. A second propagation vector,? , gives the direction of relative tangential 

1 ® 

ABAQUS is manufactured by Dassault Systemes Simulia Corp. (DSS), Providence, RI, USA 
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motion of the cohesive surfaces. At microscopic length scales where the behavior of 
discrete atoms and slip planes are significant, there is a distinction between Mode II and 
Mode III sliding and fracture. Fracture by definition generates new free surface area, 
while sliding involves bond breakage which reform between the upper and lower 
surfaces. The determination of relative amounts of tangential sliding and fracture is an 
issue that is currently unresolved and will not be elaborated here; a further discussion can 
be found in Saether (2008). 

The directional decohesion elements developed in the current work consist of a 1-D 6- 
node decohesion element for joining 2-D quadratic triangular continuum elements and a 
2-D 12-node decohesion element for joining 3-D quadratic tetrahedral elements. For both 
elements, the orientation of slip planes with respect to the cohesive surface is required 
input. To apply Rice’s criterion, the direction of Mode I opening for each decohesion 
element must be determined. 

For the 1-D decohesion element, the direction of opening or sliding is simply given 
by the sign of the slope of the relative normal displacements or relative tangential 
displacements along the element. Figure 10 shows the 6-node decohesion element 
configuration. Although the element is quadratic, the Mode I opening direction is 
determined by computing an approximate linear gradient of the relative opening 
displacements, Av„ over the element using only the end nodes. The degree of sliding in 
Mode II is computed by a summation of all relative ^-displacements, A Uj, over the 
element. 


y’> v 



Figure 10. Relative opening displacements for Mode I and Mode II fracture in 
a 6-node quadratic decohesion element. 
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The kinematics and integration scheme used in the new 6-node 1-D directional 
decohesion element is similar to the element presented by Alfano and Crisfield (2001). 
The sign of the difference between relative displacements determines the directionality of 
Mode I opening and Mode II sliding across the element. In determining the sense of 
relative motion, only the end nodes of the upper or lower surface need to be considered. 
The propagation vectors are then computed as 


d = Sgn (Av 3 - Av, )ej 

(17) 

s = Sgn (Auj-i- Au 3 )ej 

(18) 


where Sgn indicates the Sign function. 

For the 2-D 12-node decohesion element shown in Figure 11, the surface may be 
defined by the local mapped x ’ and y ’ coordinates corresponding to the node positions, 
and can be defined using a general scalar field representation as <f> = F(x \y ’), where <p 
represents the variation of a relative displacement component over this region. Using this 
surface, the components of the propagation vectors can be expressed in terms of the 
element x ’ and y ’ coordinates. The equation of a surface is given by 

Ax’ + By’ + Cz’ = D (19) 

where x’ and y’ are the mapped coordinate axes, and z’ is the value of the scalar field at 
each (x’,y’) point over the domain. The normal to this surface is given by 

n = Ae x , + Be y , + Ce z , (20) 

where e a are unit basis vectors in the local element coordinate system, and the 
components of interest are calculated as 
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= Yi (<l >2 - <J> 3 ) + Y 2 (4*3 -4>i) + y 3 (4>i -4> 2 ) 


A = 


i y'i 4>i 
i y 2 4> 2 
i y 3 4> 3 


B = 


Xi 1 f 

x 2 1 4> 2 

xi 1 (|>3 


= <I>1 (x 2 -X 3 ) + (l ) 2 (x 3 -X 1 ) + (|)3 (Xj -X 2 ). 


( 21 ) 


All global displacement quantities are mapped to a local element coordinate system for 
the determination of opening behavior. Failure propagation vectors are computed in the 
local frame and mapped back to the global coordinate system to compare with grain 
orientation angles to select which CZM law is to be applied. 

The 12-node decohesion element shown in Figure 1 1 is defined with respect to a local 
(x\y\z') element coordinate system, and the kinematics are similar to an element 
presented by Segurado and Llorca (2004). 



e i 

Figure 1 1 . Estimated direction of GB opening and sliding in 
mixed-mode fracture. 
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For Mode I, the expected opening direction is given by the propagation vector, d , as 


d 


opening 


— Aj e t + Bj e 2 


( 22 ) 


where the A f and Bi are given in Equations (21) with fa = Aw,. Because the current 
development assumes a symmetric tilt GB, the orientation angle 0 specifies the relative 
angle between a slip plane and the cohesive zone crack plane. The direction of opening in 
Mode I is given by (3 = tan -1 B ; / A r . 


This element can be enhanced to include a full description of grain boundary 
orientation with respect to surrounding slip planes. For the present analysis, it is assumed 
that the crack surface is offset from the nearest slip plane by a single tilt angle, 0 as 
shown in Figure 9. Thus, defining a local coordinate system, (x’,y’), aligned with the GB, 
the direction of opening is simply given by the sign of the Aj component of the 
propagation vector d where Sgn(d/) < 0 and Sgn(d/) > 0 indicates opening in the -x and 
+x directions, respectively. 


The component directions of Mode II and Mode III fracture are given by Equations 


(21) with fa = A Ui and fa = Av„ respectively. The propagation vectors for these two 
modes are given by 


s 


Au 


®Av 


- A n ej + B n e 2 
= A in ej + B in e 2 


(23) 


The direction of the in-plane propagation vectors may then be used to quantitatively 
determine the degree of mode mixity for fracture tangential to the crack plane. 
Directional dependence of Mode II and Mode III fracture at atomistic length scales have 
yet to be developed, however, once defined with suitable selection criteria, these may be 
included into the element formulation in an identical manner as with the automated 
selection of appropriate Mode I CZMs. 
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4.0 Simulation of Fracture in 2-D and 3-D Aluminum Polycrystals 


Investigations of fracture in 2-D and 3-D aluminum poly crystals have been 
performed. One study involves a 2-D polycrystal model in which CZM parameters 
describing brittle and ductile intergranular fracture propagation are varied to simulate the 
effect of different GB fracture behavior. This simulation does not include directionality 
of fracture but is instead focused on the effect of CZM parameterization on the overall 
stress-strain behavior of the model and final fracture path. Another study examines two 3- 
D polycrystal models to assess the effectiveness of directional decohesion finite elements 
that automatically select a brittle or ductile CZM using Rice’s criterion based on GB 
orientation with respect to adjacent slip planes and the opening direction of Mode I 
fracture. 


4. 1 Study of CZM Parameter Variation on Fracture Development 

in a 2-D Polycrystal 


Figure 12 shows a tessellated 2-D poly crystalline microstructure having many grain 
boundaries. The model has dimensions of l x = l y = 47.5 nm, and consists of 50 grains. 
Six-node triangular finite elements are used to represent the elastic continuum having the 



Figure 12. Idealized microstructural model. 


elastic properties of aluminum (£'=72 
GPa, v=0.34). Parametric studies were 
performed by linearly combining the 
CZM parameters using the MD 
extracted cohesive zone curves 
describing brittle and ductile Mode I 
fracture shown in Figure 8 (Glaessgen 
et al., 2006). The polycrystalline 
microstructure is simulated under a 
uniform tensile strain U y loaded 
incrementally to 22% strain. CZM 
properties for shear sliding under 
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Mode II loading were arbitrarily selected to simulate an elastic-perfectly plastic behavior. 
All parameters are shown in Table 2. 


Using brittle Mode I CZM properties for all the GBs, Figure 13 shows the sequence 
of decohesion processes under normal applied displacements. The CZM parameters for 
Mode I opening obtained from MD analysis yielded an initial GB stiffness that is similar 
to the stiffness of the adjoining grains, and thus does not act as a numerical penalty 
constraint. Therefore, these interfaces begin to exhibit relative displacement at load 
levels below the peak traction in the elastic region of the bilinear CZM. This is depicted 
in Figure 13a where no magnification has been applied to the relative displacements of 
the cohesive zones along the GBs and the black regions indicate significant relative 
displacements along the grain interfaces prior to full opening of the CZMs. At higher 
applied strains, the lowest energy associated with GB opening dictates the location at 
which a dominant microcrack begins to form (Figure 13b). At a critical applied strain, 
the local coalescence of GB opening forms a dominant microcrack as shown in Figure 
13c. It is also seen that due to unloading after the dominant microcrack has fractured the 
model, the stiffness in other decohesion elements acts to pull the opening cohesive zones 
back to an unloaded state. 



(a) Initial opening of GBs, 
= 1 6 % 


(b) Local crack formation, 

Sy = 17% 


(c) Dominant microcrack, 
= 22 % 


Figure 13. Dominant microcrack formation in a polycrystal model using cohesive 
zone models. 
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Prior to demonstrating the behavior of directional decohesion elements, the 
qualitative influence of CZM parameters will be illustrated through comparative 
simulations in which properties for Mode I fracture are assigned as a parametric series of 
proportions of the Mode I brittle and ductile CZM parameters presented in Table 3. These 
parameters, tj>, were obtained as 


^combined tX (^brittle (1 - Ct) (^ductile 


(24) 


where a ranged from 1.0 to 0.0. The Mode II shear parameters were arbitrarily selected 
to simulate an elastic-perfectly plastic behavior. The values chosen were f mx = 0.8 GPa 
and A c mx =100 nm which were kept constant in all simulations. Material properties for 
the grains were assumed as linear isotropic with a Young’s modulus of 100.0 GPa and 
Poisson’s ratio of 0.34. 


Table 2. CZM parameters defining Mode I behavior at the 2-D crack tip. 


CZM 

Brittle Mode I 

Ductile Mode I 

f mx (GPa) 

5.00 

3.95 

A° (nm) 

1.07 

0.85 

A c mx (nm) 

2.17 

5.65 

G (J/m 2 ) 

5.43 

11.2 


Table 3. CZM parameters used to simulate different proportions of elastic 
and plastic Mode I opening behavior at the 2-D crack tip. 


Mode I CZM parameters. 


CZM 

a 

A mx (nm) 

G (J/m 2 ) 

Cl 

1.00 

2.17 

5.43 

C2 

0.75 

3.04 

7.20 

C3 

0.50 

3.91 

8.75 

C4 

0.25 

4.78 

10.1 

C5 

0.00 

5.65 

11.2 


Figure 14 shows the resulting stress-strain response of the poly crystal model using 
variations in the proportion of brittle and ductile Mode I cohesive zone behavior. The 
global stress versus strain results demonstrate the imposed trend of decreasing maximum 
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stress with increasing ductility in the CZM law. Because the relative tangential sliding 
motion in shear is assumed to be independent of Mode I separation, all the parameterized 
models follow a similar extended softening response that continues until ultimate failure. 
Due to the different parameters used to define Mode I cohesive properties of the GBs, 
three different crack paths where exhibited by the models and are shown in Figure 14. 

2.5E-09 


2.0E-09 


1.5E-09 

£ 

z 

(/) 

</> 

Q) 

5) 1.0E-09 


5.0E-10 


0.0E+00 

0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 

Strain (%) 

Figure 14. Parametric study of Mode I stress/strain response in a polycrystal 
with brittle vs ductile GBs. 



4.2 Study of Directional Decohesion Elements on Fracture 
Development in 3-D Poly crystals 

Two 3-D polycrystal models have been developed to illustrate the influence of CZM 
parameters and the functionality of directional decohesion elements that account for grain 
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orientation and Mode I fracture opening direction. The first model focuses on the 
preferred propagation of fracture along GBs forming a junction between multiple grains. 
In a 3-D model, grain boundaries typically form triple lines and quadripole junctions. 
However, to best illustrate the directional behavior of fracture propagation, the model is 
limited to 3 grains forming a triple junction. Modeling the cohesive surfaces placed along 
each GB using the 12-node directional decohesion element discussed in Section 3.1 will, 
in general, require five independent angles to fully define the orientation of a flat GB in a 
defect-free material and the availability of accurate CZMs to describe directional fracture 
characteristics. This information will be required as input to each decohesion element 
defining a particular GB. However, because of the limited characterization of GB fracture 
behavior, for the present analysis, each GB is assumed as a symmetric tilt GB that can be 
characterized by a single misorientation angle. In order to apply Rice’s criterion, the 
local direction of opening must be computed and the angle, 0, between the opening 
direction and the local slip plane orientation must be determined. 

The grains are assumed to be isotropic and linear elastic (E = 100 GPa, v= 0.3), and 
the CZM parameters used in this simulation are presented in Table 4. The directional 
decohesion element detailed in Section 3.1 was coded into a user-defined element for use 
in ABAQUS. This element routine implemented Rice’s criterion which was 
automatically applied during the simulation to determine whether brittle or ductile Mode 
I CZM parameters were to be applied. The model geometry, loading, and results of 
changing the tilt angle for each GB is shown in Figure 15. 


Table 4. CZM parameters defining GB properties in the triple junction model. 


CZM 

Brittle Mode I 

Ductile Mode I 

Shear 

^ mx (GPa) 

5.00 

8.00 

0.80 

A° (nm) 

1.07 

1.71 

0.40 

A mx (nm) 

2.17 

8.00 

1.00 


As shown in Figure 15, simulation of a 3-D grain junction under a constant biaxial 
load state shows that different fracture paths are predicted using the formulated 
directional decohesion element by changing the tilt angle associated with each GB. 
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S = 3.0 nm 




a) = 0, e 2 - %!2, e 3 = % b) 6 1 = tc/2, 6 2 = 0, d 3 = x 

Figure 15. Simulation of a triple junction showing different opening behavior 
with different assumed GB tilt angles using directional decohesion elements. 


The second model is a 3-D poly crystal model that is used to illustrate the effect of 
changing CZM parameters associated with GBs on the global stress-strain response and 
fracture surface creation. This model was constructed by extruding the 2-D model used 
in Section 4.1 in the thickness (z-) direction and is depicted in Figure 16. The same 
loading and boundary conditions shown in Figure 12 were used, however, to fix the 
model against motion in the thickness direction, roller supports were placed on the z = 0 
surface of the model. The model incorporates directional decohesion elements along each 
GB and assigns an initial random material orientation for each grain. CZM parameters are 
systematically varied to simulate the fracture behavior of other types of symmetric tilt 
GBs. 
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Figure 16. Model of a 3-D 50-grain microstructure. 

To simplify the geometry of the model, each grain was assumed to have a principal slip 
system normal to the global x-y coordinate system. Thus, grain orientation of the //’ grain 
can be assigned using a single angle, /?„, ranging from (-tz,+tz) and defined with respect to 
the x-axis as shown in Figure 17. 



Figure 17. Geometric description of grain orientation and grain boundaries. 
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For a simple tilt grain geometry, the grain boundary orientation is defined by the angle, 
at, which is measured counterclockwise to the global y-axis. This angle is computed as 


a k = Tan 


-l 


v x -w x 


v v y 


w 


y J 


(25) 


where v and w denote vertex positions defining the ends of the GB. It follows that the tilt 
angles, $>i and 6 ^ 2 , of the GB with respect to the grain slip system are given by 


e bi = Pi - V 2 + a k 
e b2 = V 2 - Pj - a k 


(26) 


A further simplification is required to allow the parameterized MD results for a 
symmetric tilt grain boundary to be used. Therefore, each tilt angle is replaced by an 
average value, 0 , obtained as 


9 = (8 m + 8 b ,)/2 = (p, -Pj)/2 (27) 

The 3-D 25-grain microstructure model, shown in Figure 16, has x and y dimensions 
of 47.5 nm and a thickness of 4.75 nm. The grain properties were modeled as linear 
isotropic with a modulus of 100 GPa and a Poisson’s ratio of 0.34. Although similar 
geometry and grain properties where used for the 2-D and 3-D poly crystal models, CZM 
parameters were not made to coincide because these parameters were used differently in 
the two models to generate effective properties for the decohesion elements in the two 
models that precludes a direct comparison of deformation response. The Mode I CZM 

parameters used z° mx = 5 GPa for brittle fracture and z° mx = 4 GPa for ductile fracture in 

the current simulation. Both CZMs used an initial stiffness of 4.69 GPa. For Mode II and 
Mode III deformation fields, the maximum traction was held constant at 0.8 GPa, and the 
maximum opening displacement was held constant at 2.0 nm. The varied CZM 
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parameters used to represent other GB configurations with differing degrees of elastic 
and plastic response are presented in Table 5. These parameters were derived from CZM 
properties calculated by Yamakov et al. (2006) for a symmetric tilt Z99 GB in aluminum. 


Table 5. CZM parameters used to simulate different proportions of elastic and 
plastic Mode I opening behavior at the 3-D crack tip. 


Brittle Mode I CZM. Peak traction held constant at f mx = 5.0 GPa. 


CZM 

A mx (nm) 

G (J/m 2 ) 

D1 

2.17 

5.42 

D2 

2.89 

7.22 

D3 

4.33 

10.84 

D4 

6.50 

16.25 


Ductile Mode I CZM. Peak traction held constant at f mx = 4.0 GPa. 


CZM 

A mx (nm) 

G (J/m 2 ) 

D1 

3.75 

7.50 

D2 

5.00 

10.0 

D3 

6.25 

12.5 

D4 

10.0 

20.0 


Figure 18 shows a representative evolution of a dominant microcrack in the 3-D 
polycrystal model. As directional decohesion elements along various GBs exhibit 
softening due to damage accumulation, local load redistribution influences which GBs 
will ultimately fracture and create part of the dominant microcrack. For simplicity in the 
CZM formulation and finite element implementation, unloading is assumed to follow a 
path back to the origin to avoid storing and enforcing as a nodal constraint any evolving 
inelastic displacements during simulation. Thus, as various regions away from the 
dominant microcrack are unloaded, the directional decohesion elements along GBs in 
these regions pull the opening GBs back to a closed state. 
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The global stress-strain response of the model is shown in Figure 19. For each set of 
CZM parameters listed in Table 5, the initial random assignment of grain orientation 
causes the slight difference in the initial elastic stiffness between ductile and brittle 
CZMs such that each model exhibits a similar global stress-strain response until local 
softening due to damage accumulation commences at -9.35% strain. The variation in the 
maximum ultimate displacement for each model shows the expected increase in the 

degree of softening with increasing A c mx in the nonlinear portion of the global stress- 

strain curves. The difference in predicted crack paths is due to changes in the internal 
load transfer due to variations in the Mode I CZM properties of the GBs. 





s= 9.68% 


f=0.0% 


£=4.40% 





£= 10 . 1 % 


£= 11 . 1 % 


£= 11.7% 


Figure 18. Fracture sequence in the 3-D polycrystal model. 
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Figure 19. Stress-strain curves for the 3-D polycrystal model. 


5. Concluding Remarks 

The directional decohesion element developed herein adds an additional feature to 
decohesion element formulations that allows the element to alter its cohesive behavior 
during the course of nonlinear finite element analysis. For simulation of microstructural 
failure within polycrystals, the directional decohesion element incorporates Rice’s 
criterion for dislocation emission to automatically apply an appropriate brittle or ductile 
CZM. The selected CZM governs intergranular fracture based on surrounding grain 
orientation and opening direction in various ffacture/sliding modes. Two simplistic 
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illustrative examples of how results from atomistic simulations might be incorporated 
within the direction decohesion elements were presented. 

For arbitrary GB orientations, the directional decohesion element will require 
additional rotation angles to be input to define the GB and perform more involved 
calculations of the relative orientation of the cohesive plane and surrounding slip 
systems. Some additional derivations might be required to modify the opening criterion 
for the case when the crack plane does not lie along an atomic plane and intersects a slip 
plane in a skewed orientation. However, the most important requirement for the realistic 
simulation of failure in metallic polycrystals is the determination of CZMs to quantify 
GB fracture behavior for the full multi-dimensional space of GB configurational 
parameters. 
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